﻿Imports WINX32Lib 'Remember to add r3eference first (Add Reference option). For MS Visual Studio .NET 2003: before adding a reference to the project, corresponding assembly should be added to the Add Reference dialog box (http://msdn.microsoft.com/en-us/library/wkze6zky%28v=vs.71%29.aspx).
Imports System.IO
Imports System.Math


Module Grain_nalysis

    Sub Analyze_Frame_Grains(ByRef Frame_To_Analyze As Array, ByRef Frame_DimX As Integer, ByRef Frame_DimY As Integer, _
                     ByRef Frame_Index As Integer, _
                     ByVal Max_Multiplicity As Integer,
                 ByVal low As Short, ByVal high As Short, ByVal Min_Intensity As Short, ByVal low1 As Short, ByVal high1 As Short, _
                 Optional ByRef Spectrum(,) As Single = Nothing, _
                 Optional ByRef Good_Events_Sum()(,) As Single = Nothing, _
                 Optional ByRef Total_Number_Of_Good_Events(,) As Integer = Nothing, _
                 Optional ByRef Intensity_Histograme()() As Integer = Nothing, _
                 Optional ByRef Multiple_Hits()(,) As Single = Nothing, _
                 Optional ByRef Total_Number_Of_Multiple_Hits(,) As Integer = Nothing, _
                 Optional ByRef toc2 As String = "")

        Dim tic As String
        tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60

        If Not Frame_To_Analyze.Rank = 2 Then
            MsgBox("The frame to analize hase a dimention of " & Frame_To_Analyze.Rank)
        End If

        Dim Max_Intensity As Integer
        Dim Energy_Window_Center As Single = (high + low) / 2
        Dim Energy_Window_Center1 As Single = (high1 + low1) / 2
        Dim Energy_Window_HW As Single = (high - low) / 2
        Dim Energy_Window_HW1 As Single = (high1 - low1) / 2


        Dim Lower_Levels(200) As Integer
        Dim Lower_Levels1(200) As Integer

        Dim Upper_Levels(200) As Integer
        Dim Upper_Levels1(200) As Integer

        For n = 1 To 200 - 1
            Lower_Levels(n) = n * Energy_Window_Center - Sqrt(n) * Energy_Window_HW
            Lower_Levels1(n) = n * Energy_Window_Center1 - Sqrt(n) * Energy_Window_HW1
            Upper_Levels(n) = n * Energy_Window_Center + Sqrt(n) * Energy_Window_HW
            Upper_Levels1(n) = n * Energy_Window_Center1 + Sqrt(n) * Energy_Window_HW1
        Next

        Dim Frame_As_Grain_Numbers(Frame_DimX, Frame_DimX) As Integer
        Dim Hits_above_threshold(Frame_DimX, Frame_DimX) As Boolean

edges_clear_out:
        'We need to clear out the edges (alternatively we could put our array into a bigger one with zeros around)
        For i = 0 To Frame_DimX
            Frame_To_Analyze(i, 0) = 0
            Frame_To_Analyze(i, Frame_DimY) = 0
        Next
        For j = 1 To Frame_DimY - 1
            Frame_To_Analyze(0, j) = 0
            Frame_To_Analyze(Frame_DimX, 0) = 0
        Next


        ' The length of the grains array
        Dim NumEl = 10000 'Frame_DimX * Frame_DimY - 1

        Dim actual_grain_number As Integer = 1
        Dim grain_intensity() As Integer
        Dim number_of_pixels_in_grain() As Integer
        Dim number_of_hits_in_grain() As Integer
        Dim Grain_Mutiplicity(,) As Integer
        Dim Grain_Scale_Factor(,) As Double
        ReDim grain_intensity(NumEl)
        ReDim number_of_pixels_in_grain(NumEl)
        ReDim number_of_hits_in_grain(NumEl)
        ReDim Grain_Mutiplicity(NumEl, 2)
        ReDim Grain_Scale_Factor(NumEl, 2)


        'Here we check wich pixels should be analize -
        ' we choose only the ones that have the intensity abowe the treshold 
        ' and with  the sum of naighbourf abowe low1er limit of the energy window 
        For i = 1 To Frame_DimX - 1
            For j = 1 To Frame_DimY - 1
                If Frame_To_Analyze(i, j) > 0 Then
                    'Histogram is calculated here
                    Intensity_Histograme(0)(Frame_To_Analyze(i, j)) += 1
                    If Frame_To_Analyze(i, j) > Min_Intensity Then
        
                        Dim Neighbours_Sum As Integer
                        Neighbours_Sum = _
                         Frame_To_Analyze(i - 1, j - 1) + _
                         Frame_To_Analyze(i - 1, j) + _
                         Frame_To_Analyze(i - 1, j + 1) + _
                         Frame_To_Analyze(i, j - 1) + _
                         Frame_To_Analyze(i, j) + _
                         Frame_To_Analyze(i, j + 1) + _
                         Frame_To_Analyze(i + 1, j - 1) + _
                         Frame_To_Analyze(i + 1, j) + _
                         Frame_To_Analyze(i + 1, j + 1)
                        Intensity_Histograme(2)(Neighbours_Sum) += 1
                        If Neighbours_Sum > Lower_Levels1(1) Then
                            Hits_above_threshold(i, j) = True
                            Frame_As_Grain_Numbers(i - 1, j - 1) -= 1
                            Frame_As_Grain_Numbers(i - 1, j) -= 1
                            Frame_As_Grain_Numbers(i - 1, j + 1) -= 1
                            Frame_As_Grain_Numbers(i, j - 1) -= 1
                            Frame_As_Grain_Numbers(i, j) -= 2 'to distiguish from borders
                            Frame_As_Grain_Numbers(i, j + 1) -= 1
                            Frame_As_Grain_Numbers(i + 1, j - 1) -= 1
                            Frame_As_Grain_Numbers(i + 1, j) -= 1
                            Frame_As_Grain_Numbers(i + 1, j + 1) -= 1
                        ElseIf Frame_To_Analyze(i, j) > Lower_Levels(1) _
                            AndAlso Frame_To_Analyze(i, j) < Upper_Levels(1) Then
                            'just in case:
                            Frame_To_Analyze(i, j) = 0

                            'output:
                            Spectrum(i, 0) += 1
                            Good_Events_Sum(0)(i, j) += 1
                            Total_Number_Of_Good_Events(Frame_Index, 0) += 1
                            Multiple_Hits(0)(i, j) += 1
                            Total_Number_Of_Multiple_Hits(Frame_Index, 0) += 1
                        End If
                    End If
                End If
            Next
        Next


        For i = 0 To Frame_DimX
            For j = 0 To Frame_DimY
                If Hits_above_threshold(i, j) = True AndAlso Frame_As_Grain_Numbers(i, j) < 0 Then
                    '############### grain (domain) definition ##############
                    Frame_As_Grain_Numbers(i, j) = actual_grain_number
                    ' Sum of intensities:
                    Dim Sum = 0
                    Dim Good_pixels_Sum = 0
                    'idexes of pixels above the threshold:
                    Dim Good_pts As New Stack(Of Integer)(1000)
                    ' Start with the original point in the stack.
                    Dim pts As New Stack(Of Integer)(1000)
                    'We push i then j, remember to pop j then i
                    pts.Push(i)
                    pts.Push(j)
                    ' While the stack is not empty, process a pixel.
                    Do While pts.Count > 0
                        Dim jj As Integer = pts.Pop()
                        Dim ii As Integer = pts.Pop()
                        Sum += Frame_To_Analyze(ii, jj)
                        number_of_pixels_in_grain(actual_grain_number) += 1
                        ' we can stop if we are sure that we are at the border of a grain
                        If Not Frame_As_Grain_Numbers(ii, jj) = -1 Then
                            ' if we are at the hit position we check all the pixels around 
                            ' (here we do not need to check if we are not going out of Frame)
                            If Hits_above_threshold(ii, jj) Then
                                Good_pixels_Sum += Frame_To_Analyze(ii, jj)
                                Good_pts.Push(ii)
                                Good_pts.Push(jj)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii - 1, jj - 1, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii - 1, jj, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii - 1, jj + 1, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii, jj - 1, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii, jj + 1, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii + 1, jj - 1, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii + 1, jj, actual_grain_number)
                                Check_If_Pixel_Is_In_Grain(Frame_As_Grain_Numbers, pts, ii + 1, jj + 1, actual_grain_number)
                            Else
                                ' we only look at diagonal pixels 
                                ' (here we need to check if we are not going out of Frame)
                                If ii > 0 Then
                                    If jj > 0 Then _
                                        Check_If_Pixel_Is_A_Hit(Frame_As_Grain_Numbers, Hits_above_threshold, pts, ii - 1, jj - 1, actual_grain_number)
                                    If jj < Frame_DimY - 1 Then _
                                        Check_If_Pixel_Is_A_Hit(Frame_As_Grain_Numbers, Hits_above_threshold, pts, ii - 1, jj + 1, actual_grain_number)
                                End If
                                If ii < Frame_DimX - 1 Then
                                    If jj > 0 Then _
                                        Check_If_Pixel_Is_A_Hit(Frame_As_Grain_Numbers, Hits_above_threshold, pts, ii + 1, jj - 1, actual_grain_number)
                                    If jj < Frame_DimY - 1 Then _
                                        Check_If_Pixel_Is_A_Hit(Frame_As_Grain_Numbers, Hits_above_threshold, pts, ii + 1, jj + 1, actual_grain_number)
                                End If
                            End If
                        End If
                    Loop
                    grain_intensity(actual_grain_number) = Sum
                    number_of_hits_in_grain(actual_grain_number) = Good_pts.Count / 2
                    '#################### end of grain definition #########

                    Max_Intensity = Max(Max_Intensity, Sum)
                    Intensity_Histograme(1)(Sum) += 1

                    Dim Lower_Lewel_Multiplicity As Integer = Check_Lower_Limit(Sum, Energy_Window_Center, Lower_Levels1)
                    Dim Upper_Lewel_Multiplicity As Integer = Check_Upper_Limit(Sum, Energy_Window_Center, Upper_Levels1)


                    Dim multiplicity_Delta As Integer = Upper_Lewel_Multiplicity - Check_Lower_Limit(Sum, Energy_Window_Center, Lower_Levels1)

                    'here we check either we are in an energy window
                    If Lower_Lewel_Multiplicity >= Upper_Lewel_Multiplicity Then
                        If Upper_Lewel_Multiplicity = 1 Then
                            Dim max_pixel As Integer = 0
                            Dim i_max_pixel As Integer = 0
                            Dim j_max_pixel As Integer = 0
                            For k = 1 To number_of_hits_in_grain(actual_grain_number)
                                Dim jj As Integer = Good_pts.Pop()
                                Dim ii As Integer = Good_pts.Pop()
                                If Frame_To_Analyze(ii, jj) > max_pixel Then
                                    i_max_pixel = ii
                                    j_max_pixel = jj
                                    max_pixel = Frame_To_Analyze(ii, jj)
                                End If
                            Next
                            Spectrum(i_max_pixel, 0) += Upper_Lewel_Multiplicity
                            Good_Events_Sum(0)(i_max_pixel, j_max_pixel) += Upper_Lewel_Multiplicity
                            Total_Number_Of_Good_Events(Frame_Index, 0) += Upper_Lewel_Multiplicity
                            Multiple_Hits(0)(i_max_pixel, j_max_pixel) += Upper_Lewel_Multiplicity
                            Total_Number_Of_Multiple_Hits(Frame_Index, 0) += Upper_Lewel_Multiplicity
                            'If Lower_Lewel_Multiplicity = 2 Then
                            '    Spectrum(i_max_pixel, 1) += multiplicity_Delta
                            '    Good_Events_Sum(1)(i_max_pixel, j_max_pixel) += multiplicity_Delta
                            '    Total_Number_Of_Good_Events(Frame_Index, 1) += multiplicity_Delta
                            'End If
                        ElseIf Upper_Lewel_Multiplicity <= Max_Multiplicity Then
                            For k = 1 To number_of_hits_in_grain(actual_grain_number)
                                Dim jj As Integer = Good_pts.Pop()
                                Dim ii As Integer = Good_pts.Pop()
                                Spectrum(ii, 0) += Frame_To_Analyze(ii, jj) * Upper_Lewel_Multiplicity / Good_pixels_Sum
                                '     Spectrum(ii, 1) += Frame_To_Analyze(i, j) * multiplicity_Delta / Good_pixels_Sum
                                Good_Events_Sum(0)(ii, jj) += Frame_To_Analyze(ii, jj) * Upper_Lewel_Multiplicity / Good_pixels_Sum
                                '      Good_Events_Sum(1)(ii, jj) += Frame_To_Analyze(i, j) * multiplicity_Delta / Good_pixels_Sum
                                Multiple_Hits(Upper_Lewel_Multiplicity - 1)(ii, jj) += Frame_To_Analyze(ii, jj) * Upper_Lewel_Multiplicity / Good_pixels_Sum
                            Next
                            Total_Number_Of_Good_Events(Frame_Index, 0) += Upper_Lewel_Multiplicity
                            '   Total_Number_Of_Good_Events(Frame_Index, 1) += multiplicity_Delta
                            Total_Number_Of_Multiple_Hits(Frame_Index, Upper_Lewel_Multiplicity - 1) += 1
                        Else
                            For k = 1 To number_of_hits_in_grain(actual_grain_number)
                                Dim jj As Integer = Good_pts.Pop()
                                Dim ii As Integer = Good_pts.Pop()
                                Spectrum(ii, 0) += Frame_To_Analyze(ii, jj) * Upper_Lewel_Multiplicity / Good_pixels_Sum
                                Spectrum(ii, 1) += Frame_To_Analyze(ii, jj) * multiplicity_Delta / Good_pixels_Sum
                                Good_Events_Sum(0)(ii, jj) += Frame_To_Analyze(ii, jj) * Upper_Lewel_Multiplicity / Good_pixels_Sum
                                Good_Events_Sum(1)(ii, jj) += Frame_To_Analyze(ii, jj) * multiplicity_Delta / Good_pixels_Sum
                                '       Multiple_Hits(Upper_Lewel_Multiplicity - 1)(i, j) += Frame_To_Analyze(i, j) * multiplicity_Delta / Good_pixels_Sum
                            Next
                            Total_Number_Of_Good_Events(Frame_Index, 0) += Upper_Lewel_Multiplicity
                            Total_Number_Of_Good_Events(Frame_Index, 1) += multiplicity_Delta
                            '    Total_Number_Of_Multiple_Hits(Frame_Index, Upper_Lewel_Multiplicity - 1) += 1
                        End If
                    Else 'Domains with no multiplicity
                        For k = 1 To number_of_hits_in_grain(actual_grain_number)
                            Dim jj As Integer = Good_pts.Pop()
                            Dim ii As Integer = Good_pts.Pop()
                            Lower_Lewel_Multiplicity = Check_Lower_Limit(Frame_To_Analyze(ii, jj), Energy_Window_Center, Lower_Levels)
                            Upper_Lewel_Multiplicity = Check_Upper_Limit(Frame_To_Analyze(ii, jj), Energy_Window_Center, Upper_Levels)
                            If Lower_Lewel_Multiplicity = Upper_Lewel_Multiplicity Then
                                If Upper_Lewel_Multiplicity = 1 Then
                                    Spectrum(ii, 0) += 1
                                    Good_Events_Sum(0)(ii, jj) += 1
                                    Total_Number_Of_Good_Events(Frame_Index, 0) += 1
                                    Multiple_Hits(0)(ii, jj) += 1
                                    Total_Number_Of_Multiple_Hits(Frame_Index, 0) += 1
                                ElseIf Upper_Lewel_Multiplicity = 2 Then
                                    Spectrum(ii, 0) += 2
                                    Good_Events_Sum(0)(ii, jj) += 2
                                    Total_Number_Of_Good_Events(Frame_Index, 0) += 2
                                    Multiple_Hits(1)(ii, jj) += 2
                                    Total_Number_Of_Multiple_Hits(Frame_Index, 0) += 2
                                End If
                            End If
                        Next
                    End If
                    actual_grain_number += 1
                End If
            Next
        Next
        toc2 = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
    End Sub


    ' See if this point should be added to the stack.
    Private Sub Check_If_Pixel_Is_In_Grain(ByVal Frame_As_Grain_Numbers(,) As Integer, ByVal pts As Stack(Of Integer), _
        ByVal x As Integer, ByVal y As Integer, _
         ByVal new_value As Integer)
        If Frame_As_Grain_Numbers(x, y) < 0 Then
            Frame_As_Grain_Numbers(x, y) = new_value
            pts.Push(x)
            pts.Push(y)
        End If
    End Sub

    Private Sub Check_If_Pixel_Is_A_Hit(ByVal Frame_As_Grain_Numbers(,) As Integer, ByVal Hits_above_threshold(,) As Boolean, ByVal pts As Stack(Of Integer), _
    ByVal x As Integer, ByVal y As Integer, _
     ByVal new_value As Integer)
        'We do nothing if we are sure to be on the grains border (Frame_As_Grain_Numbers(x, y) = -1 ) or if we are out of the grain Frame_As_Grain_Numbers(x, y) >= 0
        If Frame_As_Grain_Numbers(x, y) < -1 Then
            If Hits_above_threshold(x, y) Then
                Frame_As_Grain_Numbers(x, y) = new_value
                pts.Push(x)
                pts.Push(y)
            End If
        End If
    End Sub


    Private Function Check_Lower_Limit(ByVal Value As Integer, ByVal Energy_Window_Center As Single, ByVal Lower_Levels As Integer()) As Integer
        Check_Lower_Limit = Value \ Energy_Window_Center
        If Value > Lower_Levels(Check_Lower_Limit + 1) Then
            Check_Lower_Limit += 1
        End If
    End Function

    Private Function Check_Upper_Limit(ByVal Value As Integer, ByVal Energy_Window_Center As Single, ByVal Upper_Levels As Integer()) As Integer
        Check_Upper_Limit = Ceiling(Value / Energy_Window_Center)
        If Value < Upper_Levels(Check_Upper_Limit - 1) Then
            Check_Upper_Limit -= 1
        End If
    End Function

End Module

